A classical theory for second-harmonic generation from metallic nanoparticles 
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In this article, we develop a classical electrodynamic theory to study the optical nonlinearities 
of metallic nanoparticles. The quasi-free electrons inside the metal are approximated as a classi- 
cal Coulomb-interacting electron gas, and their motion under the excitation of an external elec- 
. . , tromagnetic field is described by the plasma equations. This theory is further tailored to study 

0^ ' second-harmonic generation. Through detailed experiment-theory comparisons, we validate this 

' classical theory as well as the associated numerical algorithm. It is demonstrated that our theory 

' not only provides qualitative agreement with experiments, it also reproduces the overall strength of 

(N : the experimentally observed second-harmonic signals. 

PACS numbers: 42.70.-a, 52.35.Mw 
0\ ■ I. INTRODUCTION 

(N : 

I— I, Optical second-harmonic generation (SHG) from a metal (silver) surface was first observed in 1965 pLi], four years 
Q ' after the first observation of SHG from quartz in 1961 0]. In the following fifty years, a number of important features 
• i-H of SHG from metallic surfaces have been founded such as (1) second harmonic (SH) intensities can be enhanced 
more than an order of magnitude by coupling incident light into surface polariton resonances at metal surfaces 
O ! 0! (2) SHG from surfaces of centrosymmetric metals is anisotropic, the strength of the SH response thus depends 
^ . on the relative orientation of the incident field and the crystal axes [J, [^; (3) Because of local-field enhancement, 
O SHG is very sensitive to surface roughness and chemical processes such as adsorption and electrochemical reactions 
i [1, 0|- On the theoretical side, different approaches on both phenomenological and microscopic levels have been 
J>-»' developed to analyze SH response from metals 0, such as classical Boltzman equation approach [§], hydrodynamic 
'. model [13, mi 
susceptibility 
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14, 1 15 1, phenomenological formalism in terms of the fundamental tensor elements of the SH 
17, lis 



19|, |20|, |2l|, [24I and the self-consistent density functional formalism (see Ref. [23|, |2J| and 
\ the cited references). 

fSJ ■ Recently, a renaissance of scientific interests appears in the quadratic nonlinearities of metallic nanostructures 
K«*" ' and nanoparticles (NPs) partially owing to the significant localizations of electromagnetic (EM) field induced by the 
. plasmonic oscillations of the conduction electrons inside the metal [H, IH, [l^, IH, |23, [13, HH, [H, [H, HI, [H, [H, [13, 
' IS, m, EO, El, 113, m, 0, m. Hi, 1131 . More specifically, SHG were experimentally observed from different geometric 
, configurations such as sharp metal tips psi . [33 | , periodic nanostructured metal films [29j , imperfect spheres [sil [35| , 
• ■ split-ring resonators [34, .40] and their complementary counterparts ^42] , metallodielectric multilayer photonic-band- 
] gap structures [4l[ , T-shaped H^l and L-shaped NPs [H, HI] , nonccntrosymmetric T-shaped nanodimers [H, \^ and 
"fishnet" structures [ii ]. 

It should be emphasized that these subwavelength NPs and one-dimensional interfaces have different structural 
symmetries, and these differences further lead to significant consequences. For ideally infinite interfaces, the dominant 
SH dipole source appears only at the interface between centrosymmetric media where the inversion symmetry is broken, 
although higher order multipole sources provide a relatively small bulk SH polarization density. The SH polarization 
5^ ■ density is thus significantly localized in a surface region a few Angstrom wide, and sensitively influenced by the 
5^ \ details of the surface electronic structure. On the other hand, for low-symmetric or even asymmetric NPs, such as 
gold split-ring resonators, SH dipolar polarizability may be presented in the whole volume and not limited at the 
interface [i^. Consequently, the overall shape of the NP plays a significant role in determining the SH response. To 
analyze the quadratic nonlinearities of these metallic NPs, we therefore propose that complicated microscopic models 
of the interfaces are not required, and an easy-to-implement classical model is sufficient. 

The paper is organized as follows. Section 2 presents a classical electrodynamic model which describes the nonlin- 
earities induced by Coulomb-interacting electron gas in metals. Using small nonlinearity approximation, this classical 
model is further tailored in Section 3 to treat second-order generation. Section 4 gives a detailed comparison between 
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theoretical results and the corresponding experiments. Discussion, conclusion and acknowledgement are presented in 
Section 5 and Section 6, respectively. 



II. CLASSICAL ELECTRODYNAMIC MODEL 



In our model, the motion of quasi-free electrons inside a metal is described classically. Quantum mechanical 
Coulomb correlations and exchange contributions are thus missing while the classical Coulomb interaction (i.e. the 
Hartree term in a quantum mechanical derivation) is fully included [49| . Furthermore, the Coulomb scattering due to 
higher-order quantum corrections is phenomenologically described via an inverse decay time 7 = 1/r. The electrons 
inside the metal are described via their number density rip and velocity Ug. In our classical model, these quantities 
are continuous functions of position r and time t [sol, lEll, IH, HI] . We further assume the mass of ions arc infinite. 
Consequently, the ionic density no(r) is time- independent and only the electrons can move and contribute to the 
current density. In other words, an infinite barrier surface potential is assumed and no(r) is constant within the metal 
and zero outside the metal. 

We begin with two equations for electronic number density ?^e(^, t) and the velocity field Ue(r), 

^ + V • (neUe) = 0, (1) 

me + Ue ■ Ue ^ ~e (E + X B) . (2) 

Here, the first equation is the usual continuity equation expressed in terms of carrier density instead of charge density. 
The second equation is the generalization of Newton's equation to the case of a continuous field. The term in brackets 
on the left-hand side is the so-called convective or material derivative, which is a result of the description of electrons 
in terms of a continuous density and can be formally derived from the quantum mechanical Wigner distribution [4^ . 
More intuitively, it can be understood as the time derivative of an electron taken with respect to a coordinate system 
which is itself moving with velocity Ue(r(t), t), given by [5l| 



~dr ~ ~dt 



\dv{t) 1 




dt 





'd_ 

dt 



(3) 



To describe the interaction between the classical electron gas and the external EM fields self-consistently, we couple 
Eqs. (HI) and Eqs. ^ to Maxwell's equations by defining the charge density and the current density 

p{v,t) = e [no(r) - ne(r,i)] , (4) 
j(i",*) = -ene(r,t)ue(r,t) = [p(r,i) - eno(r)] Ue(r,i), (5) 

in terms of the electronic number density and velocity field. Using these definitions and the equations of motion, 
Eqs. ((D) and Eqs. ([2]), we achieve 



dp 
dt 



V-j, (6) 



, , . E [pE+jxB]-7j, (7) 

dt ^ drk \ eno - p J rrie me 

where we have added a phenomenological term —jj to describe the current decay due to Coulomb scattering. The 
Lorentz force describes a change in momentum due to an applied force. The first term on the right-hand side, resulting 
from the convective derivative, describes an increase or decrease of momentum simply due to an accumulation or 
depletion of electrons at a certain spatial point. 

Equations (O and ([7]) have to be coupled to Maxwell's equations. The final full set of equations to be solved by a 
numerical scheme read as 

^^-VxE, (8) 

I = ^E-7j + E^f-^)--[/'E+jxB], (10) 
dt rUe ^ drk \ eno - pj "le 
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where the charge density p has to be viewed as a function of the electric field since each occurrence of p can be 
replaced by the relation 

p = eoV-E. (11) 

This set of equations couples the dynamics of the EM field to the dynamics of the carriers described by their current 
density j. It should be mentioned that Equation (|10[) contains rich physics. The first two terms represent the linear 
collective oscillation of the electrons with respect to the background ionic density no(r), and the last two terms 
introduce three distinct sources for nonlinearities of the plasma. The second and third source are the well-known 
electric- and magnetic-component of the Lorentz force, respectively. The first source term is a generalized divergence 
originating from the convective time derivative of the electron velocity field Ug mentioned above. 



III. PERTURBATIVE EXPANSION OF NONLINEARITIES 



In order to obtain a simplified set of equations more suitable for a numerical approach we expand every quantity 
in a power series of the peak electric- field amplitude |i?oxc| of the excitation pulse. Formally, we can write 



E(r,t) - ^E«(r,0, 

j 

B(r,t) = ^B(j")(r,i), 

j 

j(r,t) = $]j(^)(r,t). 



(12) 
(13) 
(14) 



where the functions E^^^, B^-'^ and j^^^ scale like liJcxcP- A similar expansion automatically holds for the charge 
density by inserting Eq.(12) into Eq.(ll) 



(15) 



Separating different orders, we obtain the linear response of the metal via 

9B(i) 



dt 
9EW 

dt 
dt 



= -V X e(^\ 
= X B(i) - -J 



(1) 



= -7J 



(1) ^ ^Z!£e(i) 



(16) 
(17) 
(18) 



This is equivalent to the well-known Drude model of a metal whose bulk plasma frequency uf^ = e^no/nieeo, as can 
be easily seen by Fourier transformation [s^ . 



The second-order fields describe the lowest-order nonlinearity of the metal and are given by 

5B(2) 



dt 

dE^ _ 
dt 

OP ^ 
dt 

together with the nonlinear source term 



-V X E(2), 

C2VXB(2)-Ij(2), 



£0 



(2) + £!!^E(2) + s(2), 

m. 



dri. 



eno 



e 



eo (v-E(i')E«H-j« xBW 



(19) 
(20) 
(21) 

(22) 



where k represents x, y and z coordinate. The homogeneous part of this set of equations is identical to the first-order 
equations such that the propagation of the SH field is modified by the Drude response of the metal. The source term 
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is expressed fully in terms of the first-order fields such that the sets of Eqs. (|16p - p8)) and (|T9|) - ([22)) can be solved 
separately. 

We want to stress that no approximations have been done yet except the expansion in orders of the exciting electric 
field. All fields are real quantities and no expansion in terms of "phase factor times slowly varying envelop" (see the 
Appendix) has been done so far. In principle, these equations can be numerically solved, and a switch-off analysis 
can be further utilized to distinguish the contribution of three nonlinear sources. 

A three-dimensional finite-difference time-domain (FDTD) algorithm is employed to numerically solve the first- 
order and second-order equations separately [5j|. Yee's discretization scheme is employed so that all field variables 
are defined in a cubic grid. Electric and magnetic fields are temporally separated by half a time step, they are also 
spatially interlaced by half a grid cell. Central differences in both space and time are then applied to Eqs. (|16p - (P^ 
[55| . In addition, all the NPs studied are made of gold with Drude-type permittivity approximated as 

e{Lo) = 1.0 - , . , (23) 

with the bulk plasma frequency = 1.367 x lO^^s"^ and the phenomenological collision frequency 7 = 6.478 x lO^^s^^ 
[13, 53, [13] (please notice that we employ Wp, not no, in the simulations). In order to describe the energy conversion 
efficiency in the nonlinear-optical process, we define a normalized SH intensity, 

77=|E(2)(2^o)/E(i)(a;o)|', (24) 

to measure the strength of the far-field SH signal, where wq is the frequency of the incident fundamental-frequency 
(FF) wave. 

Our interests in the present article are limited to arrays of metallic NPs (also named as planar metamaterial 
[13, SO, Hi) with normal incidence, the computational domain is therefore arranged as follows: an array of NPs is 
placed in the middle of the space with its top and bottom surfaces positioned perpendicular to the z direction; plane 
waves propagating along the z axis are generated by a total field/scattering field technique [13]; perfect matched 
absorbing boundary conditions are applied at the top and bottom of the computational space together with periodic 
boundary conditions on other boundaries jssj l; the structure studied extends periodically in the x and y directions, 
and only single unit cell is needed in the computational space. In addition, in all the following simulations, the size 
of the spatial grid cell is fixed as 2.5 nm, and the associate time step is 4.17 attosecond. 



IV. THEORY-EXPERIMENT COMPARISON 



In this section detailed comparisons between the numerical evaluations of our theory with the experiments done 
in two independent laboratories are presented [13, [H, SO, IS^] ■ First, we consider a series of experiments reported in 
Refs . (33. l40l. [4^ ■ in which NPs with different geometrical configurations are studied (see Fig. 1). Among them the U- 
shaped NPs (split-ring resonators) are known to possess negative effective permeabilities in certain frequency regions 
and are generally referred to as magnetic metamaterials [3^ ^5^, t60|, [6l|, 113 • Strictly speaking, those "metamaterials" 
are only the first step towards a true three-dimensional bulk material. So far, most of the metamaterials are rather 
two-dimensional arrays of unit cells to study the fundamental properties of the NPs. These samples are supported 
by infinite-thickness glass (with e — 2.25) substrate coated with a thin film of indium-tin-oxide (with e — 3.8), and 
the thicknesses of the gold and indium-tin-oxide layers are 25 nm and 5 nm, respectively [13, Ho, 113 • It should be 
emphasized that the geometrical parameters of these NPs are chosen such that each structure has a resonance around 
1500-nm wavelength. 

Our theoretical results are summarized in Fig. 1. We note that the simulations qualitatively agree with the 
corresponding experimental measurements. The SH signal strength emitted from the U-shape particles with x- 
polarized FF incidence is found as 6.6 x 10^^^ from the simulation, which is quite close to the experimental result of 
2.0 X 10"^^ [40] (Please notice that the SH strength reported in Ref.[40] has been corrected in the sequent erratum). 
Our simulation thus reproduces the strength of the experimental SH signal. The following important conclusions can 
be further extracted from Fig. 1: 

(1) The polarization state of the far-field SH signal is always y-polarized (Fig. 2) for both x- or y-polarized incident 
fields (except the rectangle-shaped NPs). There thus exists a universal selection rule, that is, a mirror symmetry 
of the metallic NPs in one direction completely prohibits a polarization component of SHG in that direction. This 
symmetry dependence can be explained as follows. Within the electric dipole approximation, the far-field SH electric 
field can be related to the incident FF field such that [13, [11] 



E(2c^) = Y''^ • E(^)E(l^), 



(25) 
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FIG. 1: (Color online) Comparisons of numerical simulations and experiments for different arrays of gold nanoparticles. The 
different columns (from left to right) show shape of nanoparticles, the polarization of the incident light (indicated by the arrows), 
the linear transmission (solid) and reflection (dashed) spectra, theoretical SHG spectra (amplified 13 orders of magnitude), 
the relative strengths obtained by the theory and the corresponding experiments (inside brackets). For all the structures, the 
polarization of the generated second-harmonic waves are along the y direction. The illuminating fundamental- frequency wave 
has wavelength around 1500 nm and amplitude of 2 x 10^ (V/m). (a) The U-shape particle corresponds to the experimental 
sample shown in Fig.(lA) of Ref. 34] (also at Fig. (la) of Ref. 40]), with lattice constant = ay =305 nm. (b) The C-shape 
particle corresponds to the experimental sample shown in Fig.(lC) of Ref. [3^ (also at Fig.(lc) of Ref. [iol]). with = 567.5 
nm and ay — 590 nm. (c) The inverse-U-shape particle corresponds to the experimental sample shown in Ref. 



ay =305 nm. (d) The T-shape particle corresponds to the experimental sample shown in Fig. (2c) of Ref. 



with 
with 



Ux = 295 nm and Uy = 465 nm. (e) The rectangle-shape corresponds to the experimental sample shown in Fig.(2b) of Ref. [401] . 
Ifere and ay are the lattice constants along x and y direction, respectively. All nanoparticles dimensions shown are in 
nanometer. 



6 




FIG. 2: (Color online) The polarization state of the second harmonic emission from an array of U-shape particles illuminated 
with a x-polarized plane wave at the fundamental frequency. The second-harmonic signal is a function of the measuring angle 
(not the polarization of the incident field) 9 {0 = corresponds to the x direction) . The corresponding experiment measurement 
is shown in Fig.(2B) of Ref.[3|]. 



where *x*^^'' stands for a generalized dyadic second-order nonlinear susceptibility. The x-coordinate mirror symmetry 

(2) (2) (2) (2) 

results in the following vanishing elements, Xxxx, Xxyy, Xyxy and Xyyx- Further representing the FF field as E(z,Ci;) = 
j^^^i{ut-kz) [cos^ej: -I- sinOey], with k = uj/c being the wave number, 6 being the polarization angle and (e^) being 
the unit vector along the x (y) direction, we obtain 

E,{2lu) = 

Ey{2L0) = 

The Ex{2u}) component is simply proportional to sin 20, and therefore vanishes for 9 ~ or 9 ~ 7r/2, corresponding 
to X- or y-polarized FF incidence. On the other hand, the relationship between Ey{2uj) and 9 is nontrivial, and it 
survives for 6* = and 9 = 7r/2. 

(2) Similar to the fact that SHG at metal surfaces can be significantly enhanced by coupling incident light into surface 
polariton resonances 8] , the presence of structural plasmonic resonances also can greatly enhance SHG from metallic 
NPs. Moreover, different-order plasmonic resonances make different contributions to the SHG. For example, an 
enlarged version of the U-shape NP from Fig. (la) possesses a second-order resonance coincident with the fundamental 
resonance of the original. The SHG emitting from the fundamental resonance of the original U is considerably stronger 
than the SHG from the second-order resonance of the larger structure, even without the perfect phase matching 
requirement, because the near field enhancement is maximized for the fundamental resonance [25l [SSl. [63l. [63| . 

(3) It is found in the simulation that no SH signal is emitted from the rectangle-shaped NPs in the far field. 
As stated earlier, this is a direct consequence of the fact that an individual rectangle-shaped NP possesses mirror 
symmetries along both x and y directions. A dipole SH source is therefore forbidden and only quadrupole sources 
(and higher-order multipoles) are allowed (the retardation effects are negligible here since the thicknesses of NPs are 
much smaller than the SH wavelength, while these effects are found to excite nonlocal SH dipole for large-size gold 
nanospheres [sT]). For a periodic array of rectangle-shaped NPs with translational symmetry, the SH signal from 
each individual NP interfere destructively. Therefore, only near-field SH signal exists for the array. However, slight 
far- field SH signal is observed in the experiment. This deviation is believed to originate from the fact that the samples 
are not rigorously inversion-symmetry because of the fabrication imperfections (see the scanning electron micrograph 
shown in Ref. [40]). 

Next we study the effect of the gap on the far-field SH strength in noncentrosymmetric T-shaped gold nanodimers 
(Fig. 3). The corresponding experiment is reported in Ref. [39| . and the scanning electron micrograph images of two 
dimers are shown in Fig. 3. Obviously, these T-shaped dimers does not possess any mirror symmetry along either x 
or y direction. The gold array is 20 nm thick, with a lattice spacing of 400 nm. It is further covered with a 20-nm 
protective layer of glass and supported with an infinite-thickness glass substrate. 



E',X^^lsm20, 



E' 
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vyxx 



cos 
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sm 



(26) 
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FIG. 3: (Color online) Second harmonic generation in noncentrosymmetric nanodimers with varying gaps. The configuration 
such as YXX indicates j/-polarized second-harmonic signal with a;-polarized fundamental field. The two images show the 
unit cells as simulated, containing two structures which were digitized from the scanning electron micrograph images of the 
experimental samples (Ref.js^). The corresponding experimental results are reported in Ref. [39|] (Fig.l and Fig. 3). 



To include the relative difference between the configurations of the samples, the NPs (for all gaps) employed in our 
simulations are obtained by directly scanning the experimental samples, and the computational cell consists of two T- 
shaped dimcrs. Our numerical results arc plotted in Fig. 3 and reproduce the experimental observations qualitatively. 
More specifically, the SHG for two configurations, YXX (y-polarized SH signal with x-polarized fundamental fields) 
and YYY (y-polarized SH signal with y-polarized fundamental fields), strongly depend on the size of the gap, in a 
non-monotonically decreasing fashion. The 40-nm gap yields weak SHG responses for both configurations and the 
largest SHG response occurs for YXX from the 2-nm gap. For gaps smaller than 15 nm, the YYY response is much 
weaker than the YXX response. On the other hand in a simulation of ideal structures, without the geometrical 
variation of the dimers induced by the fabrication imperfections, the YXX response decreases monotonically with 
increasing gap. The non-monotonic responses observed in the experiments thus arise from two sources, the near-field 
enhancement around the NPs decreases with increasing the gap, and variations of the overall shape of different-gap 
nanodimers due to the fabrication imperfections. 



V. DISCUSSION AND CONCLUSIONS 



The experiment-theory agreement presented above suggests that, in contrast to an ideally infinite interface whose 
SH response strongly localizes in its surface region, the overall shape of the NP plays an important or even dominant 
role in determining efficient SH emission [48l] . Hence, even without an accurate description of the surface electrons, 
our classical model can not only successfully reproduce the experimental observations qualitatively, but also reproduce 
the SH intensities. 

It should be mentioned that Schaich developed an approach quite similar to ours to study SHG by periodically- 
structured metal surface s 165 L The major difference is that analytical parametrizations of SHG at the (top) interface 
of a thick metal slab [U l66l . 67l . [ssj are taken by assuming the parametrization scheme still works even when the flat 
metal surfaces are not of infinite extent but have edges and corners. However, the validation of this assumption is 
unclear, especially for subwavelength objects where the separation between two neighboring edges may only be tens 
of nanometers and rapid variations of the parametrizations are therefore expected. Furthermore, this parametrization 
scheme limits the application of his approach for NP with complicated boundaries such as the nanodimers studied 
above as well as NPs with thin thicknesses. In addition, there are some generalized theoretical works regarding 
nonlinear metamaterials [gl, [t^ and the nonlinear properties of negative-index metamaterials (tT], [t^ . 

We want to point out that our classical model contains only the influence of the conduction electrons and neglects 
contributions from other sources such as core electrons and lattice phonons. It therefore, for example, cannot correctly 
describe third-order nonlinearities where the electronic polarization is negligible comparing with other effects such as 
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saturated atomic absorption |52l |. To include the third-order nonlinearities, we need to add another current term 

jeW = x(')||EpE, (27) 

where x*^'^^ is the third-order nonlinear susceptibility, which equals 7.56 x 10~^^m^/V^ for gold [s^. The new set of 
equations has been utilized to study third-harmonic generation from the NPs reported in Ref. [40| . It was demonstrated 
(see Appendix II) that our simulations not only reproduce the overall strength of the experimentally observed third 
harmonic signals, but also qualitatively reproduce the structure dependent changes. As expected, we found the third 
harmonic strength to be closely related to the localization degree of the FF field inside the metal. 

In conclusion, a classical theory of second-harmonic generation from metallic nanoparticles is presented. The 
conductor-band electrons inside the metal are approximated as a classical continuous plasmonic fluid, and its dy- 
namics under an external electromagnetic field are described by the plasma wave equations self-consistently. A 
three-dimensional finite-difference time-domain approach is further applied to solve these equations numerically. By 
comparing theoretical results directly with the corresponding experiments, it is demonstrated that our classical theory, 
even without an accurate treatment of the surface electrons, qualitatively captures the dominant physical mechanisms 
of second-harmonic generations from metallic nanoparticles. This agreement suggests that the second-harmonic emis- 
sion from nanoparticles depends strongly on their overall configurations. 
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VII. APPENDIX I: APPROXIMATION FOR QUASI-MONOCHROMATIC EXCITATION 



For a quasi-monochromatic pulse with central angular frequency uq, one can classify the different contributions in 
terms of their complex phase factor. For example, the linear electric field is given by 



E(i)(r,i) = E(i)(r,t)e 



with the slowly varying complex field E'"'^^ , while the second order field 



E(2)(r,i) = E(2)(r,t)+ [Ef (r,i)e 



-i2ujQt 



(28) 



(29) 



has a second-harmonic contribution proportional to the phase factor e ^^'^o* multiplied with the slowly varying complex 

— (2) . — (2) 

amplitude Ej , as well as a slowly varying low-frequency part Eq . The magnetic field and the current can be 
expanded in a similar way. 

As a next step, we want to express the source term from Eq. ([22l) solely in terms of the linear electric field. To that 
aim, we use the linear Eqs. ()16|) and ()18p with the quasi-monochromatic approximation of Eq. p8|) and obtain 



iwoB(i) V X E(i) 



LOq 



■(1) 



'13 



'(1) ^ ^e(i) 



r(i) 



(30) 
(31) 



where we have matched the terms with equal phase factor e"*'^"*. 

Since every contribution to S^^^^ in Eq. ((22|) is of the form of a product A^^^B^^^ between two first order terms, these 
products according Eq. (|28p can be expressed as 



-iujQt 



= (1(1) exp --+C.C 
= 1(1)^(1) exp-^2'^«* 



-l-c.c. 



) exp-*"°* -l-c.c. 



A(i)(5(i))* + 



c.c. 



(32) 
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Thus, for the SH source , only the products of the slowly varying complex fields have to be calculated. They can 
be computed term by term and in the limit 7 = the first contribution fi'om the convective term is given by 



j(2) 



E 



l.k 



drk euQ 



where the plasma frequency is defined as i^pi(r) = e^no(r)/(meeo). The second term of Eq. 



(33) 



the electric Lorentz 



force - is already expressed solely in terms of the electric field and the third, magnetic term can be written as 



5(2) 



n me me UJn 2 V / 



Adding up all three contributions to the complex SH source term 83 is then given by 

,2 



o(2) _ e eo 

me UJ^ 



E 



(1) (v . (, 



2e(i) V-E(i) 



^vfE«.E(i: 



Furthermore, from the first order wave equation, we find that 

V-E(i) = ^V-(^2g(i)) 



(34) 



(35) 



(36) 



such that the SH source can be accordingly simplified to 



g(2) ^ eep 
^ me 



2EW (v.EW)+i^v(EW.E(i 



(37) 



— (2) 

In a similar fashion, also the low-frequency source Sq can be derived. Repeating analogous steps for the second 
term of Eq. ([32)1 we obtain 



g(2)^ e6o^^|g(l)|2 
me UJ„ 



(38) 



which is the well-known ponderomotive force. 

In order to insert the nonlinear source into the differential equation for j*^^), Eq. (HI]), we have to express the total 
real source in terms of the slowly varying complex amplitudes, 



S(2) 



S(2) 



g(2)g-i2a)ot 



£f£iVv|E(i)|2-i^ 

me LJq me 



2EW (v.EW)+i^v(EW.E(i 



(39) 



This result cannot be expressed by the real linear electric field for all frequencies. But since we are most interested in 
the second harmonic generation, we can approximate the source by 



S(2) 



SHG 



eep 
me 



2E(i) fv-E(i)Ui^^V|E«|2 

V / 2 



(40) 



where E*^^^ is again the full, fast oscillating, real- valued electric field obtained by the set of Eqs. (|16p ~ P^ . By inserting 
the expansion from Eq. ((28)) into Eq. (|40l) it can be easily shown that the second-harmonic contribution of Eq. ([39]) 
is exactly reproduced while the low-frequency contribution of Eq. (|^(I| is different from that of Eq. ([M)) .[73| 

To numerically solve the j2 equation with the FDTD approach, Eqs. (|19p - (PT|) with the source given by Eq. (HD|) 
have to be solved. Technically, the current is split into three different contributions according to 



Sl-(2) 2 

dt ^■'^ me 

dt 



(2)„2i^EW fv-EW) , 



(41) 
(42) 
(43) 



where the sum of j + js + jc defines the total current. 
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TABLE I: Third harmonic strength (10"^^) 



Structure 


Experimental results 


Numerical results 


U-shaped (Fig.(la)) 


3.0 


1.14' 


T-shaped (Fig. (Id)) 


0.66 


0.80 


Rectangle (Fig.(le)) 


0.21 


0.26 



"The fundamental incidence is x-polarized. 

'The difference between experiment and simulation is possibly due to experimental imperfections in the real structure, see the Fig. (2a) 



VIII. APPENDIX II: NUMERICAL RESULTS OF THIRD-HARMONIC GENERATION 

Third- harmonic generation from the samples described in Ref.fioj are numerically simulated, and the obtained 
results are listed in Table I. We see that our simulations reproduce the overall strength of the experimentally observed 
third-harmonic signals, and qualitative reproduce the structure dependent changes. As expected, we find the third- 
harmonic strength to be closely related to the localization degree of the fundamental field inside the metal. In 
addition, although the strength of second-harmonic and third-harmonic signals are comparable J3\ , almost negligible 
interactions are observed in our simulations. Second- harmonic generation and third-harmonic generation from metallic 
nanoparticles can be therefore studied separately. 
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